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Abstract. We perform high-statistics Monte Carlo simulations of three three- 
dimensional Ising spin glass models: the ± J Ising model for two values of the disorder 
parameter p, p = 1/2 andp = 0.7, and the bond-diluted ±J model for bond-occupation 
probability pb = 0.45. A finite-size scaling analysis of the quartic cumulants at the 
critical point shows conclusively that these models belong to the same universality 
class and allows us to estimate the scaling-correction exponent u> related to the leading 
irrelevant operator, lo = 1.0(1). We also determine the critical exponents v and rj. 
Taking into account the scaling corrections, we obtain v = 2.53(8) and 77 = —0.384(9). 
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The most peculiar aspect of critical phenomena is the universality of the asymptotic 
behaviour in a neighborhood of the critical point. In experiments and Monte Carlo 
(MC) simulations the possibility of approaching the critical point (and/or the infinite- 
volume limit) is generally limited. Therefore, an accurate determination of the universal 
critical behaviour requires a good control of the scaling corrections. This is particularly 
important for systems with disorder and frustration, such as spin glasses, where severe 
technical difficulties make it necessary to work with systems of relatively small size. 
Even though the critical behaviour of Ising spin glass models has been much investigated 
numerically in the last two decades, see, e.g., [HI2lSlllll5lEl[3EllSllIDlinilI2], itis 
not yet clear how reliable the numerical results are. For instance, the estimates of the 
critical exponents have significantly changed during the years, as shown by the results 
reported in pQ. Moreover, the most recent MC studies, see, e.g., [HE], find significant 
discrepancies among the estimates of the correlation-length exponent v obtained from 
the finite-size scaling (FSS) at T c of different observables, such as the temperature 
derivatives of of the Binder cumulant, and of the susceptibility. For instance, 
for the bimodal Ising model [1] quotes v = 2.39(5), v = 2.79(11), and v = 1.527(8), 
from the analysis of these three quantities. These differences indicate the presence of 
sizeable scaling corrections. However, the data are not precise enough to allow for 
scaling corrections in the analyses. In order to reduce their effects, [2] has proposed 
an alternative purely phenomenological scaling form inspired by the high-temperature 
behavior, which affects only the analytic scaling corrections and apparently reduces the 
differences among the estimates of v. Some attempts to determine the nonanalytic 
scaling corrections have been reported in [TOj [9j [8j H], but results are quite imprecise. 

Summarizing, little is known about the scaling corrections in Ising spin glass 
models, even though it is now clear that their understanding is crucial for an accurate 
determination of the critical behavior. Here we report a MC study, which represents 
a substantial progress in this direction. Indeed, by a FSS analysis of renormalisat ion- 
group (RG) invariant quantities we are able to obtain a robust estimate of the leading 
scaling-correction exponent u. This allows us to analyze the MC data for the critical 
exponents taking the scaling corrections into account. Estimates obtained from different 
observables are now in agreement. 

We consider the ±J Ising model on a cubic lattice of linear size L with periodic 
boundary conditions. The corresponding Hamiltonian is 

H — — ) J J X yO'xO' y , (1) 

where o~ x = ±1, the sum is over pairs of nearest-neighbor lattice sites, and the 
exchange interactions J xy are uncorrelated quenched random variables with probability 
distribution P(J xy ) = pS(J xy — 1) + (1 — p)S(J xy + 1). The usual bimodal Ising spin 
glass model [13], for which [J xy ] = (brackets indicate the average over the disorder 
distribution), corresponds top = 1/2. For p ^ 1/2, [J xy ] = 2p—l ^ and ferromagnetic 
(or antiferromagnetic) configurations are energetically favored. A reasonable hypothesis 
is that, along the transition line separating the paramagnetic and the spin glass phase, 
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the critical behaviour is independent of p, i.e., a nonzero value of [J xy ] is irrelevant, 
as found in mean-field models [Hj. The paramagnetic-glass transition line extends for 
p* > p > 1 — p*, where [15] p* = 0.76820(4). For 1 > p > p* the low-temperature phase 
is ferromagnetic, and the transition belongs to the randomly-dilute Ising universality 
class [16]. We also consider a bond-diluted ±J Ising model with bond-occupation 
probability pb = 0.45, and equal probability for the values ± J. 

We focus on the critical behaviour of the overlap parameter q x = a x ^a x 2 \ where 
a^ x are independent replicas with the same disorder J xy . If G(x) = [(qoq x )}, we define 
the susceptibility x = 2~2 x G{x) and the second-moment correlation length £ 

e = 1 d{0) ~ d{q) (2) 

4sin 2 (g min /2) G(q) ' 

where q = (q m ™ , 0, 0), q min = 2n/L, and G(q) is the Fourier transform of G(x). We also 
define 

R^ = ^/L, U A = -J^L, U 2 2 = -^-f — , (3) 

where = ( (XL ^ ) fe )- The quantities (jHJ) are RG invariant. We call them 
phenomenological couplings and denote them by R in the following. 

Let us first summarize some basic results concerning FSS, which allow us to 
understand the role of the analytic and nonanalytic scaling corrections. We consider 
two Ising spin glass systems coupled by an interaction h^2 x q x = h^2 x a x ^a x 2 \ in a 
finite volume of linear size L. The singular part of the corresponding disorder-averaged 
free energy density JF, which encodes the critical behavior, behaves as 

•F sing (/?, h, L) = L- d F{u h V>\u t V>\ {viL^}) (4) 

= L- d f(u h U>\u t V*) + Vul L~ d -" Uu h Ly\ Ut Lyt) + ..., 

where Uh and u t are the scaling fields associated respectively with h and t ~ T — T c 
(their RG dimensions are = (d + 2 — r])/2 and y t = 1/V), and Vi are irrelevant scaling 
fields with yi < 0. The leading nonanalytic correction-to-scaling exponent to is related 
to the RG dimension of the leading irrelevant scaling field v u = Vx, uj — —y w . The 
scaling fields are analytic functions of the system parameters — in particular, of h and 
t — and are expected not to depend on L. Note also that the size L is expected to be an 
exact scaling field for periodic boundary conditions. For a general discussion of these 
issues, see [TZl [18] and references therein. In general, u t and Uh can be expanded as 

u h = hu h {t) + 0(h 3 ), u h (t) = a h + ai t + 0(t 2 ), (5) 
ut = ctt + c 02 t 2 + c 20 h 2 + c 2 ih 2 t + 0(t 3 , h\ hH), 

where we used the fact that the free energy is symmetric under h — > —h. In the 
expansion of Uh,t around the critical point h,t = 0, the terms beyond the leading ones 
give rise to analytic scaling corrections. There are also analytic corrections due to the 
regular part of the free energy; since they scale as L v ~ 2 ~ L~ 2A , they are negligible in 
the present case. The scaling behaviour of zero-momentum thermodynamic quantities 



The critical behaviour of 3D Ising spin glass models: universality and scaling corrections 4 



can be obtained by performing appropriate derivatives of .Fsmg with respect to h. For 
instance, the overlap susceptibility x — d 2 J r /dh 2 \ h=0 behaves as 

X = L 2 ^u h (t) 2 g(u t L^) + tf-^g^utL*) + ... (6) 

The FSS of the phenomenological couplings is given by 

R0, L) = r(u t L y *) + r w (u t L y *)L-» + ... 

= R* + r'(0)ct tL yt + ... + Cu L~ w + (7) 

where R* = r (0) and c w = r w (0). In the case of Ui, this behaviour can be proved by 
taking the appropriate derivatives of T with respect to h. A similar discussion applies 
to U22 an d £/L, see Sec. 3.1 of [18] for details. The exponent v can be computed from 
the FSS of the derivative R' = dR/d(3 at /3 C , or from that of the ratio x'/Xi where 
x' = dx/d/3. At T = T c , setting t — u t — in the above-reported equations, we obtain: 



R = R* + Cl L-" + . . . , (8) 

X = cL 2 -\l + Cl L-" + ...), (9) 

R' = cL 1 / v {l + c-jT" + ...), (10) 

x > = cL 2 - T > +l ' u {\ + c x L~" + ... + ai L^ u + ...). (11) 



Note that, unlike the temperature derivative R' of an RG invariant quantity, x' a l so 
presents an L~ l l u scaling correction, due to the analytic dependence on t of the scaling 
field Uh (for this reason we call it analytic correction). Since, as we shall see, 1/v ~ 0.4 
and uj ~ 1.0, the scaling corrections in x' decay slowlier than those occurring in R'. 
This makes x'/x unsuitable for a precise determination of v and explains the significant 
discrepancies observed in pp. 

Instead of computing the various quantities at fixed Hamiltonian parameters, we 
consider the FSS keeping a phenomenological coupling R fixed at a given value Rf 
pJUdB]. This means that, for each L, we determine flf(L), such that R(/3 = (3f(L), L) = 
Rf, and then consider any quantity at (3 — (3f{L). The value Rf can be specified at will, 
as long as Rf is taken between the high- and low-temperature fixed-point values of R. 
For generic values of Rf, (3f converges to (3 C as (3f — (3 C = 0(L~ l l u ), while at Rf = R*, 
cf. (JED, (3 f - f3 c = 0{L- x l v -"). One can easily show that the FSS behaviour at fixed 
Rf = R* is given by the same general formulas derived at T c . In the case of another 
phenomenological coupling R a we have 

R a (L) = R a [(3f(L), L]^R* a + c a L~" + ..., (12) 

where R* a is universal but depends on Rf. If Rf differs from R*, x an d R' (but 
not the phenomenological couplings R) also present L~ k l v corrections with amplitudes 
proportional to (Rf — R*) k . 

In the MC simulations we employ the Metropolis algorithm, the random-exchange 
method [20], and multispin coding. We simulate the ±J Ising model at p — 0.5 for 
L=3-14,16,20,24,28, at p = 0.7 for £=3-12,14,16,20, and the bond-diluted model at 
Pb = 0.45 for L=4-12,14,16. We average over a large number N s of disorder samples: 
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Table 1. Parameters for the simulations of 'the ± J model at p — 0.5 for L > 12. We use 
the random-exchange method with Np temperatures between /3 m i n and /3 max = 0.895. 
An elementary iteration is composed by n me t Metropolis full sweeps of all Np systems 
followed by a temperature-exchange attempt of all pairs corresponding to nearby 
temperatures. The length of the each run is 48nit er iterations; the first 20n.it C r are 
discarded for equilibration. The total number of Metropolis sweeps per sample and (3 
is n-tot = 48ni ter n me t ■ The CPU time refers to a single core of a dual-core 2.4 GHz 
AMD Opteron processor. 
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^iter 
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N p 




CPU time in days 
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0.54 


308 


13 


38282 


10 


600 


288 


10 
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20 
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20 
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0.5125 
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24 
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25 


2500 


3000 


32 


0.5125 


826 


28 


285 


60 


2500 


7200 


20 


0.6575 


782 



N s « 6.4 • 10 6 up to L = 12, iV s /10 3 « 2400,2000,1500,100,46,18, respectively for 
L = 13, 14, 16, 20, 24, 28 in the case of the ± J Ising model at p = 0.5. See Table Q] for 
details. Similar statistics are collected at p = 0.7, while for the bond-diluted model 
statistics are smaller (typically, by a factor of two for the small lattices and by a factor 
of 6 for the largest ones). For each L and model we perform runs up to values of 
(3 such that R^(f3,L) is approximately 0.63, which is close to the estimates of i?| of 
PQ: 0.627(4) and 0.635(9) for an Ising model with bimodal and Gaussian distributed 
couplings, respectively. We carefully check thermalization by using the recipe outlined 
in pp. Estimates of the different observables for generic values of (3 close to (3 C are 
computed by using their second-order Taylor expansion around /3 max = 0.895. We check 
the correctness of these estimates by comparing them with those obtained by using the 
Taylor expansion around the value of j3 used in the random-exchange simulation that is 
closest to /3 max - In total, the MC simulations took approximately 30 years of CPU-time 
on a single core of a 2.4 GHz AMD Opteron processor. 

We first perform a FSS analysis at fixed = 0.63. For sufficiently large L, the 
FSS behaviour of and U22 is given by (fl~2"!) . The MC estimates of t4(L) are shown 
in Fig. [T] versus 1/L. The results for the ±J Ising model at p = 0.5 and p = 0.7 fall 
quite nicely on two straight lines approaching the same point as L — > 00, indicating that 
uj 1.0. In the case of the diluted model the approach to the large-L limit is faster: 
fits give Ua(L) = Ul + eL~ e with e w 2. This indicates that c 4 «0 [see (IT2"]) ]. According 
to the RG, this implies that the leading nonanalytic scaling correction is suppressed 
in any quantity. Thus, the approach to the critical limit should be faster, as already 
noted in [3]. We fit the data to U* + cL~ e , taking e as a free parameter. Using data for 
L > L min = 8, we obtain U* 4 = 1.514(1), U* A = 1.514(2), and JJ\ = 1.513(1) for the ±J 
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Figure 1. Phcnomenological coupling Ua{L) vs L _1 . The lines are drawn to guide 
the eye. 



model at p — 0.5 and p = 0.7, and the bond-diluted model, respectively. The fits of U22 
to U* 2 + cL~ e (L min = 6) give U* 2 = 0.1477(3), U* 2 = 0.1481(8), and U* 2 = 0.1479(6), 
respectively for the ± J Ising model at p = 0.5 and p = 0.7, and the bond-diluted model. 
These results represent a very accurate check of universality. 

The analyses of £/ 4 and U22 give also estimates of uj. The most precise ones are 
obtained from the analysis of f7 4 . In Fig. [2] we show the results for uj as obtained from fits 
of U4 to Ul + c p L~ e and of fits of the difference U±(p = 0.5; L) — U^(p = 0.7; L) to bL~ e . 
To verify the stability of the results, we have repeated the fits several times, each time 
including only the data satisfying L > L min . We estimate uj = 1.0(1). As a check, we 
verify that the ratio C22/C4 is universal [c# is the scaling-correction amplitude appearing 
in (tUP]) as predicted by standard RG arguments. Fits of U(L) — U* to cL~ e , fixing 
Ul = 1.514, U* 2 = 0.148, e = 1.0 ( L min = 12) give c 22 /c A = 0.19 and c 22 /c 4 = 0.20, 
respectively for p = 0.5 and p = 0.7, which are in good agreement. 

Equation (fl2l) holds for any chosen value of Rf. On the other hand, \ and R' do 
not present 0{L~ x l v ) corrections only if Rf = R*. Thus, before computing the critical 
exponents, we refined the estimate of i?| by performing a standard FSS analysis of 
i?£ which takes into account the scaling corrections. Fixing uj = 1.0(1), we obtained 
.R* = 0.654(7), which is slightly larger than the estimates reported in [TJ [3]. For the 
±J model at p = 0.5 we also obtained f3 c = 0.908(4). Then, in order to determine the 
critical exponent z/, we computed R'^ and U' A at fixed R^j = 0.654. In Fig. [3] we show 
results for the ±J Ising model at p — 0.5, obtained by fitting R'^ to 

lni?' = a+ -lnL + bL~ € , (13) 
v 

with e = uj = 1.0(1), as a function of L min . They are quite stable and lead to the 
estimate 



v = 2.53(6)[2], 



(14) 
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Figure 2. Estimates of the leading scaling-correction exponent u>. For each I/ m i n , 
they are obtained by fitting only the data satisfying L > £ m in- The dotted lines 
correspond to the final estimate u> = 1 .0(1). 
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Figure 3. Estimates of the exponent v from the FSS analysis of and U^, for the 
±J Ising model at p = 0.5, obtained by fitting the data to (TT5)) . The dotted lines 
correspond to the final estimate v = 2.53(8). 



where the error in brackets takes into account the uncertainty on u. Since R^j = 0.654 
is only approximately equal to _R|, we may have residual L~ x l v corrections. The 
comparison with the same analysis at fixed = 0.63 shows that their effect is negligible. 
The results from fits of to (IT3"]) . shown in Fig. El are substantially consistent. For 
example, we find v = 2.41(13) for e = 1.0 and L min = 10. For the model at p = 0.7 the 
fit of R f £ (L m i n = 8) gives v = 2.54(6). Finally we consider the bond-diluted model. If 
we use e = 2 (this is the value determined from the analysis of t7 4 and U22) we obtain 
v = 2.55(6) (I/min = 8). These results are in good agreement with the estimate (fT4"j) . 

We estimate the exponent rj by analyzing the susceptibility x a t fixed R^ = 0.654. 
We fit ln% to a+ (2 — ri)hiL + hL~ e . In the case of the ±J model at p = 0.5, fixing e = 1.0, 
we obtain 77 = -0.384(1), -0.384(4), for L min = 7,10, respectively, with xV dof £ l - 
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Our final estimate is 

v = -0.384(4) [0]{5}, (15) 

where the error in brackets is related to the uncertainty on u and that in braces gives 
the variation of the estimate as R^j varies within two error bars of i?| = 0.654(7). The 
other models give consistent, though less precise results. Finally, we have checked the 
scaling behaviour of yf, which shows L~ l l u scaling corrections for any value of R^j, 
see pip . If we take them into account, the asymptotic behaviour of x' is consistent 
with the estimates of v and r\ obtained from R' and x- F° r instance, a fit of ln%' 
to a + crlnL + ciL~ 0A + c 2 L~ l gives a = 2.78(9) (L min = 8) to be compared with 
cr = 2 — rj + \ jv = 2.78(2), obtained by using our estimates of v and r\. 

In conclusion, we have characterized the scaling corrections to the asymptotic 
critical behaviour in 3D Ising spin glass models. We have shown that the analytic t 
dependence of the scaling fields gives rise to leading L~ l l u corrections ~ 0.4) in 
the FSS of some quantities. In particular, these corrections appear in the derivative x' 
at T c . This point has been apparently overlooked in earlier FSS studies. These analytic 
corrections may also be important in other glassy systems in which v is typically large. 
We have estimated the leading nonanalytic scaling-correction exponent uj from the FSS 
of the quartic cumulants, obtaining uj = 1.0(1). Finally, we have used these results 
to obtain accurate estimates of the critical exponents v and rj. An analysis of the 
MC data that takes into account the leading scaling corrections gives v = 2.53(8) and 
n = —0.384(9). Results obtained by using different observables and different models are 
consistent. This confirms that the bimodal Ising model belongs to a unique universality 
class, for any p in the range 1 — p* < p < p*, irrespective of bond dilution. 
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